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Abstract 

The performance of many-body perturbation theory for calculating ground- 
state properties is investigated. We present fully numerical results for the 
electron gas in three and two dimensions in the framework of the GW ap- 
proximation. The overall agreement with very accurate Monte Carlo data is 
excellent, even for those ranges of densities for which the GW approach is 
often supposed to be unsuitable. The latter seems to be due to the fulfilment 
of general conservation rules. These results open further prospects for accu- 
rate calculations of ground-state properties circumventing the limitations of 

standard density functional theory. 
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I. INTRODUCTION 

Many-body perturbation theory (MBPT), particularly in Hedin's GW approximation!, 
has been used extensively to calculate quasiparticle (QP) energies and spectra of a wide vari- 
ety of electron systems!. The GW method offers a simple way to determine the one-electron 
Green's function, G, from which the QP properties can be easily extracted. However, G 
also contains information about ground-state properties: the expectation value of any one- 
particle operator can be expressed in terms of G, and by using the Galitskii-Migdal formula! 
the total energy may also be obtained. Nonetheless, the capability of a Green's function 
MBPT to provide reliable ground-state energies has not been fully explored so far. The few 
available results are restricted to the spin-unpolarized homogeneous electron gas (HEG) in 
the range of metallic densities (i.e., r s = 2 ~ 5 a.u.)EHl. These investigations suggest that the 
GW approach could produce accurate electron total energies, but a deeper study is needed 
to provide an overall assessment of this issue. This is an important question because many 
of the limitations of the usual implementations of density functional theory (DFT)B, can 
be circumvented by using MBPT total-energy calculations. Indeed, well-recognized failures 
of the DFT in its usual approximations (for instance when studying van der Waals forcesi, 
chemical reactions^, defects in semiconductors^, or quasi-two dimensional systems^) are 
mainly due to the limited account of non-local effects that, on the contrary, are included 
in the GW approximation. On the other hand, quantum Monte Carlo (QMC) methods are 
being applied to more and more systems^, but they require a large computational effort. 
In this context, MBPT has to be regarded as a good candidate to supersede standard DFT 
schemes [like the local density (LDA) and the generalized gradient approximations] without 
implying a prohibitive computational task. To provide insights into the above points, in 
this Paper we present GW results for the ground-state properties of the three-dimensional 
HEG (covering a broad range of densities in both spin-unpolarized and fully spin-polarized 
phases) and of the two-dimensional HEG. To do so we have used the space-time numerical 
procedure developed by Rojas e t alB. This method permits efficient and stable eomputa- 
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tion of the full self-energy operator E and the corresponding Green's function G0, with the 
precision needed for a converged evaluation of the total energy. 



II. THEORY 

In MBPT, the Green's function and the self-energy of a system of iV electrons under a 
external potential f ex t (r) are linked through the Dyson equation 

G- 1 (u) = G 1 (u) - [S (u) + (Av - Afj) l] , (1) 

where the usual matrix operations are implied. Go (u) is the Green's function of a fictitious 
system of N non- interacting electrons under the potential v (r) + A/i, Av = fH + ^ext — (^h 
being the exact Hartree potential), and ApL is a constant that aligns the chemical potential 
of the fictitious system with the actual one, In the GW framework, X is approximated 
by 

Z(l,2) = iG(l,2 + )W(l,2) , (2) 



where the labels 1,2 symbolize space-time coordinates. W is the screened Coulomb potential 
which is related to the bare Coulomb interaction w and the polarizability P by 

W (u) = w + w P (u) W (u) . (3) 

Finally, under the GW approach we haveB 

P(l,2) = -2zG(l,2)G(2,l + ) . (4) 



Eq. [T|4J may be solved iteratively to self-consistency. We note that the choice of the fictitious 
non-interacting system in ([I]) is arbitrary because the differences arising from different GqS 
are cancelled out by the terms Av and Afi. Also, by including the shift A// we guarantee 
that at any step of the iteration, G verifies several exact properties that have to be verified by 
any realistic Green's function^, ensuring a smooth and stable convergence of the iterative 
process. Furthermore, the Hartree potential has to be updated after each iteration, and 



this is done by calculating the electron density n (r) from G. To evaluate A/i we need to 
obtain the chemical potential that is, by definition, the energy of the highest occupied QP 
energy, calculated at each iteration in terms of the self-energy by solving the QP Schrodinger 
equation. 

The above set of equations defines a conserving approximation in the Baym-Kadanoff 
senseEJ. One consequence is that the total number of particles given by the self-consistent 
GW Green's function does not change when an external perturbation acts on the system. 
Besides, it gives the right number of particles for the HEG§. The correctness of the num- 
ber of particles for an arbitrary inhomogeneous system can thus be inferred by regarding 
the system as the result of an adiabatic transform of the HEG. Another characteristic of a 
conserving approximation is the absence of ambiguities among different expressions for cal- 
culating the total energy. This suggests that a MBPT evaluation of ground-state properties 
should employ a conserving approximation such as self-consistent GW. 

However, routine GW calculations are mainly concerned with the QP properties of real 
materials and do not attempt self-consistency. Indeed, self-consistency implies a worsening 
in the description of the QP spectrum rather than an improvements^. Hence, the usual 
non self-consistent (and non-conserving) GW approach (that we shall denoted as GqWo, 
whereas GW will stand for the fully self-consistent solution of Hedin's equations) is clearly 
preferred in a spectral context. This failing of GW can be understood in terms of the spectral 
properties of G, but our present interest is very different: issues as those described in the 
previous paragraph are by far more important than the concrete shape of G. Of course these 
two aspects are not independent, but the development and application of a conserving theory 
giving at the same time an accurate description of QP spectra is an unsolved formidable 
challenge. 

In GqWo, the self-energy is approximated by Eq. || supposing that G = Go, whereas 
the screened Coulomb potential was obtained from ([|) and (|) once (i.e., W has been 
calculated in a RPA fashion). Eventually, the Dyson equation is solved taking into account 
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that Av ~ vlda, where ^lda is the LDA exchange-correlation (XC) potential. Although the 
term A^x is often neglected at this stage, in this Paper we will keep this contribution for 
the reasons explained above. Partial self-consistency (denoted as GWq) may be achieved by 
keeping the screened Coulomb potential Wo obtained from a Go Wo procedure and, hence, by 
solving only (|I|) and (|2[) iterativelyi'0. Although the GWq is not a conserving approximation, 
it gives the right number of particles for the HEG (so meeting an important physical point 
of the GW scheme)!. On the other hand, its description of QP properties seems to be only 
marginally worse than G wM- Finally, there are many other schemes for implementation of 
the GW equations, but essentially they are focused on the choice of Go for non self-consistent 
calculations, and so are meaningless in an homogeneous system. 

As mentioned above, we will apply the space-time methodll! to solve the GW equations. 
Each one is written in the most favorable spatial representation (reciprocal or real), go- 
ing from one to other using Fourier transforms. However the most important issue is the 
evaluation of the dynamical dependence on imaginary time or frequency domain. To cal- 
culate ground-state properties, a contour deformation avoids the need to obtain G for real 
frequencies. Concretely, the expectation value of any one-particle operator b is given by 
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-Im / duTr 

IT JC 



bG(u) 



(5) 



where the frequency is measured from the chemical potential fi, and C is the integral path in 
the complex frequency plane equal to the circular arc 7 from u = -ootow= — ioo together 
with the negative imaginary axis. In the same notation, the Galitskii-Migdal formula for 
the total energy reads 

1 



E = -Im / du Tr 

7T JC 



(u + h )G(uj\, (6) 



ho being the one-electron hamiltonian with potential v = v ext —^- To deal with the evaluation 
of (|5p and (|^) we write the Green's function as G = Gx + SG, Gx being the solution of the 
Dyson's equation ([I]), but substituting the full self-energy E (u>) by its frequency-independent 
part Ex (1,2) = iG (1, 2 + ) w (1, 2). Hence, the frequency integrals can be split up in two 



parts. The contribution due to Gx is evaluated analytically, whereas for the remainder the 
only non-zero contribution arises from the imaginary axis, which is amenable for numerical 
calculation. We have used Gauss- Legendre (GL) grids for imaginary times and frequencies, 
and the contributions due to points outside the GL grids are treated in an analytical fashion 
in accordance with the overall numerical procedure given in Reffi Usually, GL grids with 
128 points suffice for well-converged results (better than 1 mHa). For the HEG, matrix 
inversions are not needed, and a fully self-consistent resolution of the GW equations only 
takes typically a few seconds on a standard workstation. 



III. RESULTS AND DISCUSSION 

Using different GW schemes, we have obtained the XC energy per particle, exc (defined 
as the difference between the energies of the interacting and non-interacting systems) for 
the spin-unpolarized (( = 0) HEG (Table p] and Fig. [I]). We first compare our numerical 
results with the two of von Barth and HolmH. A small discrepancy (1 ~ 2 mHa) appears, 
but it is consistent with the error bar (about 3%) of the semi-analytical procedure carried 
out by these authorsil. Focusing on the accuracy of the MBPT procedure, we can see that 
the agreement between the essentially exact diffusion Monte Carlo (QMC)^'" and the self- 
consistent GW is almost perfect in the limit of high densities. This is not a surprise because 
the exchange is treated exactly by the GW and it is dominant in this range of densities. 
However, the quality of the GW energies is striking for intermediate and low densities, 
where the many-body effects not included in the GW framework might be expected to be 
evident. Partial self-consistency {GWq) yields slightly inferior results, though the differences 
are no more than a few mHa. The worst results (but, in any case, with errors no greater 
than 10 mHa for metallic densities) are provided by the non-self-consistent GqWq procedure. 
GqWq underestimates the total energy, and by achieving partial self-consistency, the spectral 
weight in the Green's function is blue-shifted, so increasing the total energy. After full self- 
consistency, such shift is slightly larger, but the kinetic energy is smaller than in GWq. The 



6 



presence of these two opposite trends explains the small differences between GW and GWq. 
We have also included (for comparison) the corresponding RPA values!!. Note that the RPA 
dielectric function is the same than the GqWq one, and the huge discrepancies between them 
arise from the different ways in which the evaluation of the total energy is performed. 

All the above trends also apply for the fully spin-polarized (( = 1) HEG (see Fig. 
but in this case the errors in the GW energies are marginally greater. (Although not at all 
the objective of this Paper, it is interesting to note that using the self-consistent GW, the 
paramagnetic phase is more stable than the ferromagnetic one up to r s = 15 ~ 20, in fair 
agreement with the QMC value@ of r s = 25 ~ 30, despite the low density of the system). 

Our results for the two-dimensional (2D) HEG are also shown. In 2D systems, correlation 
effects are much more important; in other words, the diagrams that are neglected in the GW 
scheme play a relevant role in these low- dimensional problems. As a consequence we cannot 
expect here extremely accurate results using the GW approximation. However, as we can see 
in Fig. 0, the GW gives energies near the QMC values@, resolving partially the inaccuracy 
of the GqWq approach and greatly improving the RPA energies. We note that in the limit 
of low densities, GWq fits the QMC results slightly better than GW. 

Finally, it is now very well known that GqWq does not give the right number of particles 
for an inhomogeneous system0I3. However there were certain doubts whether it recovers 
the right density of the HEG or not. The use of the space-time method allows us to affirm 
that GqWq does not give exactly the number of particles also in the homogeneous limit. 
Whereas the exact density and the GqWq one are indistinguishable up to r s = 4, for r s = 5 
the GqWq overestimates the density by 0.3%. The deviation increases when going into the 
low density region, being 1.7% for r s = 10, and 6.1% for r s = 20. 

In summary, we have studied the performance of the GW approximation for the evalua- 
tion of ground-state properties. The accuracy of the results is correlated with the fulfillment 
of conservation rules (that can be achieved by using a self-consistent GW scheme), and the 
approximations inherent in the GW scheme have much less importance than when calcu- 
lating QP properties. The dynamical dependences in the GW equations are easy to handle 
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using a representation in imaginary time and frequency, that may be straightforwardly gen- 
eralized to arbitrary inhomogeneous systems. Hence, the results presented here can be the 
point of departure for future accurate evaluations of ground-state properties of electron 
systems without the limitations of DFT and the complexity of QMC. 
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TABLES 
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A 


K 
J 


i n 


9D 




U.OloU 


n 97/19 
U.Z / 


n 1 /ifi/i 

U. 14:04: 


nil Q7 




U.U0414 




0.5127 


0.2713 




0.1201 




0.0344 


GW 


0.5160(2) 


0.2727(5) 


0.1450(5) 


0.1185(5) 


0.0620(9) 


0.032(1) 






0.2741 


0.1465 








GW 


0.5218(1) 


0.2736(1) 


0.1428(1) 


0.1158(1) 


0.0605(4) 


0.030(1) 


G W 


0.5272(1) 


0.2821(1) 


0.1523(1) 


0.1247(1) 


0.0665(2) 


0.0363(5) 


RPA 


0.5370 


0.2909 


0.1613 


0.1340 


0.0764 


0.0543 


-ex 


0.4582 


0.2291 


0.1145 


0.0916 


0.0458 


0.0229 


TABLE I. Minus XC 
3D homogeneous electron 


energies per particle (in Hartrees) for the spin-unpolarized phase of the 
gas obtained through several GW schemes. The second row in the GW 



entry corresponds to Ref. [6]. Also shown are the RPA results, and the QMC values from Ref. 
[23] (first row) and from Ref. [24] (second row). Parenthesis indicates the numerical uncertainty 
in the last significant figure. For reference, the exchange energy per particle, ex, is included. 
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FIGURES 

FIG. 1. XC energy per particle, £xc> f° r the spin-unpolarized 3D homogeneous electron 
gas. The essentially exact Monte Carlo results (symbols) are compared with several GW schemes 
(lines). The excellent performance of the self-consistent GW and (to a lesser extent) the partially 
self-consistent GWq, is evident. Note that the differences between several Monte-Carlo results are 
less than the symbols size. 

FIG. 2. As Fig. |l|, for the fully spin-polarized (x = 1) phase of the 3D electron gas. 

FIG. 3. As in Fig. [l], for the spin-unpolarized 2D electron gas. 
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